Given the follow set of sample data, compare the different schemes for Geometric Brownian Motion (GBM)
$$ \begin{align} \text{Today’s stock price, } S_0 &= 100 \\[5pt] \text{Time, } T &= 1 \text{ year} \\[5pt] \text{volatility, } \sigma &= 20\% \\[5pt] \text{constant risk-free interest rate, } r &= 5\% \\[5pt] \end{align} $$The following schemes will be considered for comparison
Begin by looking at the closed form solution.
The Euler-Maruyama scheme is an approximation by taking a Taylor Series expansion on the exact solution
Let,
$$
x = ( r - \frac{1}{2} \sigma^2 ) \; \delta t \;+\; \sigma \phi \sqrt{\delta t} \\
$$
Then Taylor Series expansion on the exponential with $a=0$,
$$ \begin{align} e^x &\approx e^a \left( 1 + (x-a) + \frac{1}{2} (x-a)^2 + \frac{1}{6} (x-a)^3 \;+\; ... \right) \\ &\approx 1 + x + \frac{1}{2} x^2 + \frac{1}{6} x^3 \;+\; ... \\ \end{align}\\ $$Then dropping the higher order terms $O(\delta t)$,
$$ S_{t + \delta t} \;\sim\; S_t \; ( 1 \;+\; ( \; r - \frac{1}{2} \sigma^2 ) \; \delta t \;+\; \sigma \phi \sqrt{\delta t} \; ) \\ $$Results in an approximation with an accuracy of the order $O(\sqrt{\delta t})$
The Milstein Scheme, is an inclusion of one higher order term from the Taylor Expansion done in the Euler-Maruyama scheme.
$$ exp \left\{ \; ( r - \frac{1}{2} \sigma^2 ) \; \delta t \;+\; \sigma \phi \sqrt{\delta t} \right\} \; \sim \; 1 \;+\; ( r - \frac{1}{2} \sigma^2 ) \; \delta t \;+\; \sigma \phi \sqrt{\delta t} + \underbrace{ \frac{1}{2} \sigma^2 \phi^2 \delta t }_{ \text{Milstein correction}} + O(\delta t^2) \\ $$Similarly, dropping the higher order terms $O(\delta t^2)$,
$$ S_{t + \delta t} \;\sim\; S_t \; \left( 1 \;+\; (\; r - \frac{1}{2} \sigma^2 ) \; \delta t \;+\; \sigma \phi \sqrt{\delta t} \;+\; \frac{1}{2} \sigma^2 \phi^2 \delta t \; \right) \\ \;\sim\; S_t \; \left( 1 \;+\; r \; \delta t \;+\; \sigma \phi \sqrt{\delta t} \;+\; \frac{1}{2} \sigma^2 ( \phi^2 -1) \; \delta t \; \right) \\ $$Results in an approximation with an improved accuracy of the order $O(\delta t)$.
To run the simulation, make the following assumptions:
import numpy as np
import pandas as pd
import sympy as sp
sp.init_printing()
from functools import partial
from copy import deepcopy
from time import perf_counter
from IPython.display import display_html
import plotly.express as px
import plotly.graph_objects as go
import plotly.io as pio
pio.renderers.default = 'notebook' #_connected'
from plotly.offline import init_notebook_mode
init_notebook_mode(connected=True)
from watermark import watermark
%load_ext watermark
%watermark -u -d -v -m -iv
Last updated: 2023-04-27 Python implementation: CPython Python version : 3.9.16 IPython version : 8.12.0 Compiler : GCC 7.5.0 OS : Linux Release : 5.4.0-147-generic Machine : x86_64 Processor : x86_64 CPU cores : 8 Architecture: 64bit sympy : 1.11.1 numpy : 1.24.2 pandas: 2.0.0 plotly: 5.14.1
Define the known constants
# define constants
S0 = 100
rate = 0.05
sigma = 0.2
T = 1
t = 252
delta_t = T/t
runs = 100_000
Initialize arrays for stock price and values for the random variables for repeatability across different schemes.
# initialize empty stock price array
S = np.zeros((t, runs))
S[0] = S0
# initialize random values
np.random.seed(20230419)
arr_rand = np.random.standard_normal(t*runs).reshape(t,runs)
Create a function for each scheme.
def exact(s_prev: float, r: float, vol: float, dt: float, randvar: int) -> np.array:
return s_prev * np.exp( (r - 0.5 * vol**2) * dt + vol * np.sqrt(dt) * randvar)
def euler_maruyama(s_prev: float, r: float, vol: float, dt: float, randvar: int) -> np.array:
return s_prev * ( 1 + r * dt + vol * np.sqrt(dt) * randvar)
def milstein(s_prev: float, r: float, vol: float, dt: float, randvar: int) -> np.array:
return s_prev * ( 1 + r * dt + vol * np.sqrt(dt) * randvar + 0.5 * vol**2 * (randvar**2 - 1) * np.sqrt(dt))
ls_scheme = {'exact': partial(exact, r=rate, vol=sigma),
'euler_maruyama': partial(euler_maruyama, r=rate, vol=sigma),
'milstein': partial(milstein, r=rate, vol=sigma) }
Create a wrapper function to run Monte Carlo simulation for each scheme using the same random array.
def simulate_path(approx_func, dt: float, arr_price: np.array, arr_random: np.array) -> np.array:
arr_out = deepcopy(arr_price)
for i in range(0, arr_out.shape[0]-1):
arr_out[i+1] = approx_func( s_prev=arr_out[i], dt=dt, randvar=arr_random[i] )
return arr_out
# helper function for long table display
def style_table(df: pd.DataFrame, columns: int=3, spacing: int=20):
space = "\xa0" * spacing
k, m = divmod(df.shape[0], columns)
output = space
for j in range(columns):
df_table = df[j*k+min(j, m):(j+1)*k+min(j+1, m)]
df_styler = df_table.style.set_table_attributes("style='display:inline'")
output += df_styler._repr_html_() + space
return output
Measure the time performance of the each scheme and combine into a dataframe.
ls_out = ['']*3
%%timeit -n3
ls_out[0] = simulate_path(ls_scheme['exact'], delta_t, S, arr_rand)[-1]
337 ms ± 3.23 ms per loop (mean ± std. dev. of 7 runs, 3 loops each)
%%timeit -n3
ls_out[1] = simulate_path(ls_scheme['euler_maruyama'], delta_t, S, arr_rand)[-1]
84.1 ms ± 3.43 ms per loop (mean ± std. dev. of 7 runs, 3 loops each)
%%timeit -n3
ls_out[2] = simulate_path(ls_scheme['milstein'], delta_t, S, arr_rand)[-1]
163 ms ± 12.7 ms per loop (mean ± std. dev. of 7 runs, 3 loops each)
From the observed times, it is evident that the Euler-Maruyama scheme is the fastest and is around 4 times faster than the exact solution. While the Milstein scheme took about 50% longer than Euler-Maruyama, it is still a third of the time for the exact solution.
df_sims = pd.DataFrame( ls_out, index=['exact', 'euler_maruyama','milstein']).T
df_sims.describe()
| exact | euler_maruyama | milstein | |
|---|---|---|---|
| count | 100000.000000 | 100000.000000 | 100000.000000 |
| mean | 104.983647 | 104.984039 | 104.969628 |
| std | 21.154043 | 21.149681 | 21.348863 |
| min | 41.090237 | 41.023505 | 42.191759 |
| 25% | 89.962204 | 89.954568 | 89.803186 |
| 50% | 102.846399 | 102.840868 | 102.722896 |
| 75% | 117.775404 | 117.771684 | 117.738598 |
| max | 273.659534 | 272.468921 | 290.989523 |
The different schemes appear to have similar averages but the Milstein scheme has a slightly higher standard deviation and a higher max value.
Notice the average return is around 5% which is what is expected under risk neutral conditions since the replicating portfolio must perform the same as the riskfree asset for the no arbitrage condition to hold. Also, the median is below the mean, therefore expecting right-tailed/right-skewed distributions.
Examine the histogram of the output.
ls_hist = []
for k in range(df_sims.shape[1]):
ls_hist.append( np.histogram( df_sims.iloc[:,k], bins=100))
fig1 = go.Figure()
fig1.add_trace(go.Scatter(x=ls_hist[0][1], y=ls_hist[0][0],
mode='lines', name='Exact', line=dict(shape='hv'),
hovertemplate='Exact - bin: $%{x:.3f} \t count: %{y:,}<extra></extra>'))
fig1.add_trace(go.Scatter(x=ls_hist[1][1], y=ls_hist[1][0],
mode='lines', name='E-M', line=dict(shape='hv'),
hovertemplate='E-M - bin: $%{x:.3f} \t count: %{y:,}<extra></extra>'))
fig1.add_trace(go.Scatter(x=ls_hist[2][1], y=ls_hist[2][0],
mode='lines', name='Milstein', line=dict(shape='hv'),
hovertemplate='Milstein - bin: $%{x:.3f} \t count: %{y:,}<extra></extra>'))
fig1.update_yaxes( showspikes=True) #, range=[104, 105.3] )
fig1.update_xaxes( showspikes=True, tickformat=',')
fig1.update_traces( line_width=3, opacity=0.55)
fig1.update_layout( title='Figure 1 - Histogram of estimated Stock Price', hovermode="x",
xaxis_title="Iterations",
yaxis_title="Stock price ($)",
legend_title="Scheme",
width=850, height=550)
fig1.show()
In general, all schemes are right tailed. The output of the Milstein scheme seems to have a longer tail and higher peak compared to the rest. The Euler-Maruyama scheme appears to be tracking the exact solution quite well.
Iterate over columns to calculate mean across the number of iterations using expanding mean.
for column in df_sims:
df_sims[ column + '_mean'] = df_sims[column].expanding().mean()
Look at convergence by taking the expanding mean of dataframe.
fig2 = go.Figure()
fig2.add_trace(go.Scatter(x=df_sims.index, y=df_sims.exact_mean,
mode='lines', name='Exact',
hovertemplate='Exact: $%{y:.3f}<extra></extra>'))
fig2.add_trace(go.Scatter(x=df_sims.index, y=df_sims.euler_maruyama_mean,
mode='lines', name='E-M',
hovertemplate='E-M: $%{y:.3f}<extra></extra>'))
fig2.add_trace(go.Scatter(x=df_sims.index, y=df_sims.milstein_mean,
mode='lines', name='Milstein',
hovertemplate='Milstein: $%{y:.3f}<extra></extra>'))
fig2.update_yaxes( showspikes=True, range=[104, 105.15], tickformat='.1f' )
fig2.update_xaxes( showspikes=True, tickformat=',')
fig2.update_traces( line_width=1, opacity=0.75)
fig2.update_layout( title='Figure 2 - Price vs Iterations', hovermode="x",
xaxis_title="Iterations",
yaxis_title="Stock price ($)",
legend_title="Scheme",
width=850, height=550)
fig2.show()
In this case, the difference is very small. However, Euler-Maruyama performs better than the Milstein scheme.
Compare the error of Euler-Maruyama and Milstein to the exact solution.
Define error as the difference between the scheme and the exact, $$ \varepsilon = \hat{S}_{scheme} - S_{exact} $$
df_sims['em_error'] = df_sims.euler_maruyama_mean - df_sims.exact_mean
df_sims['milstein_error'] = df_sims.milstein_mean - df_sims.exact_mean
df_sims['em_error_percentage'] = 100 * (df_sims.em_error / df_sims.exact_mean)
df_sims['milstein_error_percentage'] = 100 * (df_sims.milstein_error / df_sims.exact_mean)
fig3 = go.Figure()
fig3.add_trace(go.Scatter(x=df_sims.index, y=df_sims.em_error_percentage,
mode='lines', name='E-M',
hovertemplate='<b>E-M error (%):</b> %{y:.4%}<extra></extra>'))
fig3.add_trace(go.Scatter(x=df_sims.index, y=df_sims.milstein_error_percentage,
mode='lines', name='Milstein',
hovertemplate='<b>Milstein error (%):</b> %{y:.4%}<extra></extra>'))
fig3.update_yaxes( showspikes=True, range=[-0.12, 0.05], tickformat='.1%', dtick=0.02 )
fig3.update_xaxes( showspikes=True)
fig3.update_traces( line_width=2.5, opacity=0.85)
fig3.update_layout( title='Figure 3 - Error comparison', hovermode="x",
xaxis_title="Iterations",
yaxis_title="Error (%)",
legend_title="Scheme",
width=850, height=550)
fig3.show()
The error is significantly smaller for Euler-Maruyama while Milstein scheme is still quite good below 2%. However, it is interesting to note that the Milstein scheme does not truly converge to the exact solution. The E-M scheme is within 0.1% after around 7,000 iterations (out of 100,000 iterations).
Based on this, it appears that for this case E-M scheme offers both the speed and accuracy over the Milstein scheme.
Since the accuracy of the approximation schemes are a function of $\delta t $, consider the effect of varying $\delta t $. This is done by increasing the timesteps (partition size) from daily to $n$ multiple days which gives round (integer) numbers for a year. Increasing step size results in lower number of partitions as the time horizon is maintained to be 1 year.
$$ \text{timestep, } t = \frac{252}{n} \\ $$For $n=1, \; \delta t = \frac{1}{252}$
For $n=2, \; \delta t = \frac{1}{252 \; / \; 2} = \frac{2}{252}$
and so forth for up to $n=42$ days (two months).
The exact list is given in the table below.
days = 42
n = [ k for k in range(1,days+1) if 252 % k == 0]
# init dataframe to store results
df_deltaT = pd.DataFrame( np.zeros((len(n),5)),
columns=['n', 't','dt','em_error_perc','milstein_error_perc'])
df_deltaT['n'] = n
df_deltaT['t'] = (252/df_deltaT.n).astype('int')
df_deltaT['dt'] = T / df_deltaT.t
df_deltaT.iloc[:,:3]
| n | t | dt | |
|---|---|---|---|
| 0 | 1 | 252 | 0.003968 |
| 1 | 2 | 126 | 0.007937 |
| 2 | 3 | 84 | 0.011905 |
| 3 | 4 | 63 | 0.015873 |
| 4 | 6 | 42 | 0.023810 |
| 5 | 7 | 36 | 0.027778 |
| 6 | 9 | 28 | 0.035714 |
| 7 | 12 | 21 | 0.047619 |
| 8 | 14 | 18 | 0.055556 |
| 9 | 18 | 14 | 0.071429 |
| 10 | 21 | 12 | 0.083333 |
| 11 | 28 | 9 | 0.111111 |
| 12 | 36 | 7 | 0.142857 |
| 13 | 42 | 6 | 0.166667 |
# init empty stock price array - init to largest required size, then slice as necessary for larger timesteps
S_dt = np.zeros(( max(df_deltaT.t), runs))
S_dt[0] = S0
# initialize random values
np.random.seed(20230419)
arr_rand_dt = np.random.standard_normal(max(df_deltaT.t)*runs).reshape( max(df_deltaT.t),runs)
for k in df_deltaT.index:
df_dt = pd.DataFrame({
'exact': simulate_path( ls_scheme['exact'], df_deltaT.dt[k], S_dt[:df_deltaT.t[k]],
arr_rand_dt[:df_deltaT.t[k]])[-1],
'euler_maruyama': simulate_path( ls_scheme['euler_maruyama'], df_deltaT.dt[k],
S_dt[:df_deltaT.t[k]], arr_rand_dt[:df_deltaT.t[k]])[-1],
'milstein': simulate_path( ls_scheme['milstein'], df_deltaT.dt[k], S_dt[:df_deltaT.t[k]],
arr_rand_dt[:df_deltaT.t[k]])[-1],
})
df_deltaT.loc[k,'em_error_perc'] = (100 * (np.mean(df_dt.euler_maruyama) - np.mean(df_dt.exact))
/ np.mean(df_dt.exact)).item()
df_deltaT.loc[k,'milstein_error_perc'] = (100 * (np.mean(df_dt.milstein) - np.mean(df_dt.exact))
/ np.mean(df_dt.exact)).item()
display_html( style_table( df_deltaT, 2), raw=True)
| n | t | dt | em_error_perc | milstein_error_perc | |
|---|---|---|---|---|---|
| 0 | 1 | 252 | 0.003968 | 0.000373 | -0.013354 |
| 1 | 2 | 126 | 0.007937 | 0.000455 | -0.015436 |
| 2 | 3 | 84 | 0.011905 | -0.000245 | -0.011239 |
| 3 | 4 | 63 | 0.015873 | -0.000210 | -0.013360 |
| 4 | 6 | 42 | 0.023810 | -0.002310 | -0.005157 |
| 5 | 7 | 36 | 0.027778 | -0.002745 | -0.005314 |
| 6 | 9 | 28 | 0.035714 | -0.002499 | -0.011095 |
| n | t | dt | em_error_perc | milstein_error_perc | |
|---|---|---|---|---|---|
| 7 | 12 | 21 | 0.047619 | -0.004517 | -0.008735 |
| 8 | 14 | 18 | 0.055556 | -0.006080 | -0.007073 |
| 9 | 18 | 14 | 0.071429 | -0.008485 | -0.006405 |
| 10 | 21 | 12 | 0.083333 | -0.010361 | -0.006508 |
| 11 | 28 | 9 | 0.111111 | -0.012498 | -0.011254 |
| 12 | 36 | 7 | 0.142857 | -0.017226 | -0.010987 |
| 13 | 42 | 6 | 0.166667 | -0.017955 | -0.015620 |
Table above showing the varying $\delta t$ and error of schemes.
fig4 = go.Figure()
fig4.add_trace(go.Scatter(x=df_deltaT.dt, y=df_deltaT.em_error_perc,
mode='lines+markers', name='E-M',
hovertemplate='<b>dt:</b> %{x:.4f}<br>' + f'<b>E-M error (%):</b> ' +
'%{y:.4%}<extra></extra>'
))
fig4.add_trace(go.Scatter(x=df_deltaT.dt, y=df_deltaT.milstein_error_perc,
mode='lines+markers', name='Milstein',
hovertemplate='<b>dt:</b> %{x:.4f}<br>' + f'<b>Milstein error (%):</b> ' +
'%{y:.4%}<extra></extra>'
))
fig4.update_yaxes( showspikes=True, tickformat='.1%', dtick=0.002 ) #range=[-0.55, 0.015],
fig4.update_xaxes( showspikes=True, tickformat='0.2f')
fig4.update_layout( title='Figure 4 - Error (%) vs Iterations' +
'<br><sup>for varying dt at 100,000 iterations</sup>',
hovermode="x", xaxis_title="dt", yaxis_title="Error (%)",
legend_title="Scheme", width=850, height=550)
fig4.show()
It is expected that as $\delta t$ increases, the error would increase. The error for Euler-Maruyama scheme appears to grow somewhat linearly. While the error of the Milstein scheme seems to be improving initially but as the step size increases, the error starts to increase as well. This error turn around is most likely driven by the correction term which has a $\sigma^2$. More importantly, these errors are at the second/third decimal places onwards which is already very small. In most cases, any decision will most likely be driven by the computing resources to run these calculations.
In conclusion, from all the comparisons above, it appears that the preferred scheme will be between the exact solution and the Euler-Maruyama scheme depending on speed, low complexity and low error. For small computations, the exact solution will still be performant, however as the simulation size grows there will be a need to improve performance.
For an Asian option, the expected value of the discounted payoff under the risk neutral density $\mathbb{Q}$ is given by,
$$ V (S, t) = e^{-r(T-t)} \; \mathbb{E^\mathbb{Q}}[ Payoff] $$For call and put options, consider the following:
Using the previously provided data, $$ \begin{align} \text{Today’s stock price, } S_0 &= 100 \\[5pt] \text{Time, } T &= 1 \text{ year} \\[5pt] \text{volatility, } \sigma &= 20\% \\[5pt] \text{constant risk-free interest rate, } r &= 5\% \\[5pt] \end{align} $$
and
$$ \begin{align} \text{Strike, } E &= 100 \\[5pt] \text{Time to expiry, } (T-t) &= 1 \text{ year} \\[5pt] \end{align} $$It is important to understand that the stock price is a Martingale under the $\mathbb{Q}$ measure and the mean or drift rate is zero. This is similar to what was considered in section 1.1.
For a floating strike, the strike price is not pre-determined but instead is calculated as a function of the stock price which means it will move as the stock price moves. Clearly, this makes the option heavily path dependent as the path the stock price takes will determine the payoff.
The calculation of the strike price may take almost any form such as weighted, exponential weighted, volume weighted and so forth. In this case this analysis will compare between geometric and arithmetic means.
In general, the payoff function for an option has the form,
$$ \begin{align} P_{call} &= max(S - E, \; 0) \\[4pt] \end{align} $$and a put option,
$$ \begin{align} P_{put} &= max(E - S, \; 0) \end{align} $$ $$ \text{where } \quad S=\text{Stock price at expiry} \quad \text{and} \quad E = \text{strike price} \\ $$These are the exact payoff function for a fixed strike call and put option (European options). None of these are path dependent as the payoff only depends on the price of the underlying asset at expiry.
For options with a floating strike, whether arithmetic or geometric, continuous or discrete payoff, it will involve modifying the strike price, $E$.
First, denote the floating strike as the function $I(S)$, hence the payoff function of the form,
$$ P_{call, floating} = max(S - I, 0) \\[5pt] $$and a put option,
$$ P_{put, floating} = max(I - S, 0) \\[5pt] $$and look at the different forms of $I(S)$.
Floating arithmetic payoff $$ I_{arithmetic}(S_t) \;=\; \frac{1}{T} \int_{0}^{T} S(\tau) \; d\tau \;=\; \frac{1}{N} \sum_{i=0}^{N} S_i \\[5pt] $$
Floating geometric payoff $$ I_{geometric}(S) = \left( \prod_{i = 0}^{N} S_i \right)^\frac{1}{N} \;=\; \sqrt[N]{S_0 * S_1*...*S_N} \;=\; exp \left\{ \frac{1}{N} \sum_{i=0}^{N} \; ln \; S_i \right\} \\[5pt] $$
For continuously sampled, $N$ is the number of days to expiry whereas for the discrete case, $N$ is the number of dates sampled.
While in theory it is possible to consider the stock price in continuous time, in the real world, this is not possible. So the stock price is discretized by sampling and this may take the form of continuous sampling at some constant rate (such as daily) or discrete sampling where the stock price is sampled on specific dates or larger than daily constant interval (such as weekly or monthly).
Of course, each of these affects the strike price in different ways and more importantly, maps how the path of the underlying asset will affect the strike price and in turn, the payoff at expiry.
In the following analysis,
Use the previous defined functions and variables, begin by defining the array of underlying prices, the discount term and a function to calculate value of the option.
# Monte Carlo Price
arr_underlying = simulate_path(ls_scheme['exact'], delta_t, S, arr_rand)
strike = 100.0
tau = 1
discount_term = np.exp( -rate*tau )
dict_discrete_sampling = {'weekly': 5, 'fortnightly': 10, 'monthly': 21}
# store results for plotting
df_call = pd.DataFrame({'type':[], 'payoff': [], 'value': []})
df_put = pd.DataFrame({'type':[], 'payoff': [], 'value': []})
dict_payoff = { 'call': {}, 'put': {}, 'hist_call': {}, 'hist_put': {} }
Define functions to perform strike calculations, expected payoff and expected value of discounted payoff.
def fixed_strike(strike, size) -> np.array:
return np.repeat(strike, size)
def arith_cont(arr_price: np.array) -> np.array:
return np.cumsum(arr_price, axis=0) / np.arange(1, arr_price.shape[0]+1)[:,np.newaxis]
def geom_cont(arr_price: np.array) -> np.array:
return np.exp(( np.cumsum( np.log(arr_price), axis=0) /
np.arange(1, arr_price.shape[0]+1)[:,np.newaxis])[-1])
def arith_discrete(arr_price: np.array, n: int) -> np.array:
return ( np.cumsum(arr_price[::n], axis=0) /
np.arange(1, arr_price[::n].shape[0]+1)[:,np.newaxis])
def geom_discrete(arr_price: np.array, n: int) -> np.array:
return np.exp(( np.cumsum( np.log(arr_price[::n]), axis=0) /
np.arange(1, arr_price[::n].shape[0]+1)[:,np.newaxis])[-1])
# expected payoff and expected value of discounted payoff
def value_option( arr_price: np.array, arr_strike: float, discounting: float, option_type: str="" ) -> np.array:
if option_type.lower() == 'call':
arr_payoff = arr_price[-1] - arr_strike
elif option_type.lower() == 'put':
arr_payoff = arr_strike - arr_price[-1]
arr_payoff_filt = np.ma.filled(np.ma.masked_less_equal(arr_payoff, 0), fill_value=0)
option_payoff = np.mean( arr_payoff_filt)
option_value = discounting * option_payoff
return option_payoff, option_value
# calculate all strikes
def calc_strikes( strike, arr_price: np.array, discrete_sample: dict) -> dict:
dict_strikes = { 'fixed': fixed_strike(strike, arr_underlying.shape[1]),
'arithmetic continuous': arith_cont(arr_underlying)[-1],
'geom continuous': geom_cont(arr_underlying) }
for key,value in dict_discrete_sampling.items():
dict_strikes[f'arithmetic discrete {key}'] = arith_discrete(arr_underlying, value)[-1]
dict_strikes[f'geom discrete {key}'] = geom_discrete(arr_underlying, value)
return dict_strikes
dict_strike = calc_strikes( strike, arr_underlying, dict_discrete_sampling)
# calculate expected payoff and expected value of discounted payoff
for i,key in enumerate(dict_strike.keys()):
df_call.loc[i,:] = (key,) + value_option( arr_underlying, dict_strike[key], discount_term, 'call')
df_put.loc[i,:] = (key,) + value_option( arr_underlying, dict_strike[key], discount_term, 'put')
dfcall_styler = df_call.sort_values('value', ascending=False
).style.set_table_attributes("style='display:inline'"
).set_caption('Call Options for various strike')
dfput_styler = df_put.sort_values('value', ascending=False
).style.set_table_attributes("style='display:inline'"
).set_caption('Put Options for various strike')
display_html( "\xa0"*20 + dfcall_styler._repr_html_() + "\xa0" * 30 + dfput_styler._repr_html_(), raw=True)
| type | payoff | value | |
|---|---|---|---|
| 0 | fixed | 10.876316 | 10.345871 |
| 8 | geom discrete monthly | 6.707436 | 6.380310 |
| 7 | arithmetic discrete monthly | 6.492180 | 6.175552 |
| 4 | geom discrete weekly | 6.327508 | 6.018911 |
| 2 | geom continuous | 6.319354 | 6.011156 |
| 6 | geom discrete fortnightly | 6.308866 | 6.001179 |
| 3 | arithmetic discrete weekly | 6.106878 | 5.809042 |
| 1 | arithmetic continuous | 6.101603 | 5.804024 |
| 5 | arithmetic discrete fortnightly | 6.083495 | 5.786799 |
| type | payoff | value | |
|---|---|---|---|
| 0 | fixed | 5.892668 | 5.605279 |
| 7 | arithmetic discrete monthly | 3.762504 | 3.579004 |
| 8 | geom discrete monthly | 3.630305 | 3.453253 |
| 1 | arithmetic continuous | 3.571504 | 3.397319 |
| 3 | arithmetic discrete weekly | 3.567416 | 3.393431 |
| 5 | arithmetic discrete fortnightly | 3.544841 | 3.371957 |
| 2 | geom continuous | 3.438694 | 3.270987 |
| 4 | geom discrete weekly | 3.433377 | 3.265929 |
| 6 | geom discrete fortnightly | 3.408852 | 3.242600 |
Table above displays the expected payoff and expected value of the discounted payoff for both call and put options for the different types of strikes ordered by descending value.
fig5 = go.Figure()
for columns in df_call.T:
fig5.add_trace(go.Scatter(x=[0,1], y=df_call.T.iloc[:0:-1,columns], mode='lines+markers',
name=f'{df_call.T.iloc[0,columns]}', meta=f'{df_call.T.iloc[0,columns]}',
hovertemplate='<b>%{meta}</b> - $%{y:.3f}<extra></extra>',
line=dict(width=1.3, dash='solid'),
marker=dict(size=8, symbol='x')
))
fig5.update_yaxes( showspikes=True)
fig5.update_xaxes( showspikes=True)
fig5.update_layout( title='Figure 5 - Monte Carlo <b>Call</b> Option Estimation' +
'<br><sup>100,000 iterations - various strike types</sup>',
hovermode="x", xaxis_title="Time", yaxis_title="Payoff ($)",
legend_title="Strike type", width=850, height=475)
fig5.show()
Among the Pythagorean means, the arithmetic mean is the largest followed by the geometric mean. Therefore, it is expected that the options with arithmetic averaging will have a lower payoff (higher strike) compared to geometric averaging. As expected, the fixed strike option will have the highest value. The limited data for discrete arithmetic mean results in a lower strike and in turn higher payoff as compared to the continuous arithmetic mean. When comparing discrete sampling regimes, monthly sampling has higher payoff versus fornightly and weekly. Essentially, the less data used for sampling results in lower strike and higher payoff. In effect allowing volatility to drive prices higher than the sampled average. This line of thinking could easily be extended to include the fixed strike as well where it has no sampling.
For the fixed strike, there is no path dependence, therefore the strike is the level for which the underlying either exceeds or does not reach at expiry. Where as for the rest, the path at the sample points dictates the strike.
Comparisons which results in higher call option value in order of significance:
fig7 = go.Figure()
for columns in df_put.T:
fig7.add_trace(go.Scatter(x=[0,1], y=df_put.T.iloc[:0:-1,columns], mode='lines+markers',
name=f'{df_put.T.iloc[0,columns]}', meta=f'{df_put.T.iloc[0,columns]}',
hovertemplate='<b>%{meta}</b> - $%{y:.3f}<extra></extra>',
line=dict( width=1.3, dash='solid'),
marker=dict(size=8, symbol='x')
))
fig7.update_yaxes( showspikes=True)
fig7.update_xaxes( showspikes=True)
fig7.update_layout( title='Figure 7 - Monte Carlo <b>Put</b> Option Estimation' +
'<br><sup>100,000 iterations - various strike types</sup>',
hovermode="x", xaxis_title="Time", yaxis_title="Payoff ($)",
legend_title="Strike type", width=870, height=475)
fig7.show()
The chart above shows the expected payoff and the expected value of the discounted payoff at $t=0$ and $t=T=1$. Similar to the chart for call options, the relationship is defined by the discounting factor.
The first thing that is noticeable is that the values for the floating strike are approximately half of the fixed strike. The fixed strike has the highest value.
In the case of the put option, it is desirable to have strike as large as possible since the payoff function is, $P = I - S$. Therefore for puts, arithmetic mean results in higher value than geometric averaging.
For sampling time, it is similar to call options where less sampling results in higher value.
Comparisons which results in higher put option value in order of significance:
Previously, the effect of varying $\delta t$ has been investigated. This section will vary the data to discover the effect on option price. Specifically the volatility and time to expiry. The expected effects are:
ls_sigma = np.linspace(0.02,0.4,20) #np.linspace(0.1,0.6,11)
colnames_vol = ['vols'] + list(dict_strike.keys())
dict_vols = { 'call': pd.DataFrame( columns=colnames_vol),
'put': pd.DataFrame( columns=colnames_vol) }
dict_vols['call'].loc[:,'vols'] = dict_vols['put'].loc[:,'vols'] = ls_sigma
ls_sigma
array([0.02, 0.04, 0.06, 0.08, 0.1 , 0.12, 0.14, 0.16, 0.18, 0.2 , 0.22,
0.24, 0.26, 0.28, 0.3 , 0.32, 0.34, 0.36, 0.38, 0.4 ])
for j,vol in enumerate(ls_sigma):
arr_underlying1 = simulate_path( partial(exact, r=rate, vol=vol), delta_t, S, arr_rand)
dict_strike = calc_strikes( strike, arr_underlying1, dict_discrete_sampling)
for i,key in enumerate(dict_strike.keys()):
dict_vols['call'].loc[j, key] = value_option( arr_underlying1, dict_strike[key],
discount_term, 'call')[1]
dict_vols['put'].loc[j, key] = value_option( arr_underlying1, dict_strike[key],
discount_term, 'put')[1]
fig8 = go.Figure()
for column in dict_vols['call'].drop('vols', axis=1):
fig8.add_trace(go.Scatter(x=dict_vols['call'].vols, y=dict_vols['call'].loc[:,column], mode='lines+markers',
name=f'{column}', meta=f'{column}',
hovertemplate='%{meta}: $%{y:.3f}<extra></extra>',
line=dict( width=1.3, dash='solid'),
marker=dict(size=8, symbol='x')
))
fig8.update_yaxes( showspikes=True)
fig8.update_xaxes( showspikes=True, dtick=0.04, tickformat='.2f')
fig8.update_layout( title='Figure 8 - Call Option Value for Varying volatility' +
'<br><sup>100,000 iterations - various strike types</sup>',
hovermode="x", xaxis_title="volatility", yaxis_title="Value ($)",
legend_title="Strike type", width=870, height=600)
fig8.show()
In general, it is expected that as volatility increases, the option value increases as well. Since low volatility indicates that the underlying is less likely to move whereas a higher volatility means that the underlying will be move higher or lower at each step. Interestingly, for the floating strike options there seems to be some second order relationship at volatilities below 0.2 and beyond 0.2, the relationship tends to be linear. This behaved is observed in the fixed strike as well although much less pronounced. Potentially this is due to the riskfree term $r - \frac{1}{2} \sigma^2$ in the underlying and as $\sigma$ becomes larger, the diffusion term takes over. The floating strike in a way amplifies this as the averaging is able to keep up with volatility is low.
fig9 = go.Figure()
for column in dict_vols['put'].drop('vols', axis=1):
fig9.add_trace(go.Scatter(x=dict_vols['put'].vols, y=dict_vols['put'].loc[:,column], mode='lines+markers',
name=f'{column}',
meta=f'{column}',
hovertemplate='%{meta}: $%{y:.3f}<extra></extra>',
line=dict( width=1.3, dash='solid'),
marker=dict(size=8, symbol='x')
))
fig9.update_yaxes( showspikes=True)
fig9.update_xaxes( showspikes=True, dtick=0.04, tickformat='.2f')
fig9.update_layout( title='Figure 9 - Put Option Value for Varying volatility' +
'<br><sup>100,000 iterations - various strike types</sup>',
hovermode="x", xaxis_title="volatility", yaxis_title="Value ($)",
legend_title="Strike type", width=870, height=600)
fig9.show()
Similar to the call option, increasing volatility increases the put option value. However, it is interesting to note that the option value is almost zero as the underlying price will not have move significantly. Yet, for the Asian option floating strike, the value appears to be higher for low volatility and has the same second order hockey stick profile.
Investigate the effect of varying time to expiry while keeping all other data constant, specifically $\delta t$ shall remain at 1 day. The time to expiry being considered as in increments of 2 months from 1 months up till 18 months.
months = [ i for i in range(1, 19, 1)]
colnames_mte = ['months', 'horizon', 'step', 'discount'] + list(dict_strike.keys())
dict_mte = { 'call': pd.DataFrame( columns=colnames_mte),
'put': pd.DataFrame( columns=colnames_mte)}
for key in dict_vols.keys():
dict_mte[key].loc[:,'months'] = months
dict_mte[key].loc[:,'horizon'] = dict_mte[key].months / 12
dict_mte[key].loc[:,'step'] = (dict_mte[key].horizon * 252).astype('int')
dict_mte[key].loc[:,'discount'] = np.exp( (-rate*dict_mte[key].horizon).to_list())
dict_mte['call'].iloc[:,:4]
| months | horizon | step | discount | |
|---|---|---|---|---|
| 0 | 1 | 0.083333 | 21 | 0.995842 |
| 1 | 2 | 0.166667 | 42 | 0.991701 |
| 2 | 3 | 0.25 | 63 | 0.987578 |
| 3 | 4 | 0.333333 | 84 | 0.983471 |
| 4 | 5 | 0.416667 | 105 | 0.979382 |
| 5 | 6 | 0.5 | 126 | 0.97531 |
| 6 | 7 | 0.583333 | 147 | 0.971255 |
| 7 | 8 | 0.666667 | 168 | 0.967216 |
| 8 | 9 | 0.75 | 189 | 0.963194 |
| 9 | 10 | 0.833333 | 210 | 0.959189 |
| 10 | 11 | 0.916667 | 231 | 0.955201 |
| 11 | 12 | 1.0 | 252 | 0.951229 |
| 12 | 13 | 1.083333 | 273 | 0.947274 |
| 13 | 14 | 1.166667 | 294 | 0.943335 |
| 14 | 15 | 1.25 | 315 | 0.939413 |
| 15 | 16 | 1.333333 | 336 | 0.935507 |
| 16 | 17 | 1.416667 | 357 | 0.931617 |
| 17 | 18 | 1.5 | 378 | 0.927743 |
# initialize arrays to largest required size, then slice as necessary for larger timesteps
# initialize empty stock price array
S_mte = np.zeros(( max(dict_mte['call'].step), runs))
S_mte[0] = S0
# initialize random values
np.random.seed(20230419)
arr_rand_mte = np.random.standard_normal(max(dict_mte['call'].step)*runs
).reshape( max(dict_mte['call'].step),runs)
for j in dict_mte['call'].index:
arr_underlying2 = simulate_path( partial(exact, r=rate, vol=sigma), delta_t,
S_mte[:dict_mte['call'].step[j]],
arr_rand_mte[:dict_mte['call'].step[j]])
dict_strike = calc_strikes( strike, arr_underlying2, dict_discrete_sampling)
for i,key in enumerate(dict_strike.keys()):
dict_mte['call'].loc[j, key] = value_option( arr_underlying2, dict_strike[key],
dict_mte['call'].discount[j], 'call')[1]
dict_mte['put'].loc[j, key] = value_option( arr_underlying2, dict_strike[key],
dict_mte['call'].discount[j], 'put')[1]
fig10 = go.Figure()
for column in dict_mte['call'].iloc[:,4:]:
fig10.add_trace(go.Scatter(x=dict_mte['call'].months, y=dict_mte['call'].loc[:,column],
mode='lines+markers', name=f'{column}', meta=f'{column}',
hovertemplate='%{meta}: $%{y:.3f}<extra></extra>',
line=dict( width=1.3, dash='solid'),
marker=dict(size=8, symbol='x')
))
fig10.update_yaxes( showspikes=True)
fig10.update_xaxes( showspikes=True, dtick = 3)
fig10.update_layout( title='Figure 10 - Call Option Value for Varying Time to Expiry' +
'<br><sup>100,000 iterations - various strike types</sup>',
hovermode="x", xaxis_title="months", yaxis_title="Value ($)",
legend_title="Strike type", width=870, height=600)
fig10.show()
As expected, time to expiry is proportional to the value of the option. Given that the underlying has sufficiently time to diffuse or move through its price range. Option value is also influenced by the discounting factor relationship $e^{-rT}$ as the time to expiry increases. In the chart above, the hockey stick profile is also visible
fig11 = go.Figure()
for column in dict_mte['put'].iloc[:,4:]:
fig11.add_trace(go.Scatter(x=dict_mte['put'].months, y=dict_mte['put'].loc[:,column],
mode='lines+markers', name=f'{column}', meta=f'{column}',
hovertemplate='%{meta}: $%{y:.3f}<extra></extra>',
line=dict( width=1.3, dash='solid'),
marker=dict(size=8, symbol='x')
))
fig11.update_yaxes( showspikes=True)
fig11.update_xaxes( showspikes=True, dtick = 3)
fig11.update_layout( title='Figure 11 - Put Option Value for Varying Time to Expiry' +
'<br><sup>100,000 iterations - various strike types</sup>',
hovermode="x", xaxis_title="months", yaxis_title="Value ($)",
legend_title="Strike type", width=870, height=600)
fig11.show()
The chart above for put options indicates that for short expiry options (2 months and below) for this set of data, the floating strike scheme has high option value with high sampling arithmetic averaging being highest. Again, this is potentially due to the low data samples for averaging.
The use of closed form solution, Euler-Maruyama or Milstein scheme will depend on the use case. However, for running Monte Carlo simulation of stock prices, the Milstein scheme does not seem to offer any advantage as the Euler-Maruyama appears to offer the trifecta of speed, better accuracy and faster convergence.
The comparisons above highlight the different methods for strike calculation and it is clear in all cases that a fixed strike (European option), with all else being equal, should be preferred as there are no path dependence. However, as demonstrated, path dependence can be taken into account in the simulation.
In order of highest value for options:
Further investigation could be done for floating price options as well as other method for floating strike calculations such as volume weighted average or exponential weighted average.
The model problem is given as
$$ \frac{d^2y}{dx^2} + P\;(x) \frac{dy}{dx} + Q\,(x)\;y = f\,(x) \\[15pt] \text{subject to the boundary conditions,} \quad y\;(a) = \alpha \quad \text{and} \quad y\;(b) = \beta \\ $$and $x$ is given as a regular partition of the interval $[a,b]$ where,
$$ x_0 < x_1 < x_2 <....< x_{n-1} < x_n \\[5pt] \text{such that} \\ x_0 = a \quad \text{and} \quad x_n = b \\ $$with the partition size,
$$ \delta x = \frac{b-a}{n} $$Note that the model has a single input variable, $x$.
Let, $$ y_i \;=\; y\;(x_i) \\ P_i \;=\; P\;(x_i) \\ Q_i \;=\; Q\;(x_i) \\ f_i \;=\; f\;(x_i) \\ $$
Using central difference approximation, start by applying Taylor Series Expansion on $y\,(x)$, for $ x + \delta x $ around $ x $
$$ y\;(x + \delta x) \;=\; y\;(x) \;+\; \delta x \frac{dy}{dx} \;+\; \frac{\delta x^2}{2} \frac{d^2y}{dx^2} \;+\; \frac{\delta x^3}{3!} \frac{ d^3 y}{ dx^3} \;+\; O\;(\delta x^4) \tag{1} \\ $$and similarly for $ x - \delta x $ around $ x $
$$ y\;(x - \delta x) \;=\; y\;(x) \; - \; \delta x \frac{dy}{dx} \;+\; \frac{\delta x^2}{2} \frac{d^2y}{dx^2} \;-\; \frac{\delta x^3}{3!} \frac{ d^3 y}{ dx^3} \;+\; O\;(\delta x^4) \tag{2} \\ $$Consider the difference of $\;(1) \;-\; (2)$,
$$ \begin{align} y\;(x+ \delta x) \;-\; y\;(x-\delta x) \;&=\; 2\; \delta x\; \frac{dy}{dx} + O\;(\delta x^3) \\[8pt] \frac{dy}{dx} \;&=\; \frac{y\;(x+ \delta x) \;-\; y\;(x-\delta x)}{ 2\; \delta x\; } \\[8pt] \;&=\; \frac{ y_{i+1} \;-\; y_{i-1} }{ 2 \; \delta x } \qquad for \quad i = 1,2,...,n-1\\[8pt] \end{align} $$Consider the sum of $\;(1) \;+\; (2)$,
$$ \begin{align} y\;(x+ \delta x) \;+\; y\;(x-\delta x) \;&=\; 2y\;(x) \;+\; \delta x^2 \; \frac{d^2 y}{dx^2} \;+\; O\;(\delta x^4) \\[8pt] \frac{ y\;(x+ \delta x) \;-\; 2y\;(x) \;+\; y\;(x-\delta x) } { \delta x^2 } \;&=\; \frac{ d^2 y }{ dx^2 } \;+\; O\;(\delta x^4) \\[8pt] \frac{d^2 y}{dx^2} \;&=\; \frac{ y\;(x+ \delta x) \;-\; 2y\;(x) \;+\; y\;(x-\delta x) }{ \delta x^2 } \\[8pt] \;&=\; \frac{ y_{i+1} \;-\; 2y_i \;+\; y_{i-1} }{ \delta x^2 } \qquad for \quad i = 1,2,...,n-1\\[8pt] \end{align} $$Next, substitute the numerical approximations of the derivatives found above into the model.
$$ \frac{ y_{i+1} \;-\; 2y_i \;+\; y_{i-1} }{ \delta x^2 } \;+\; P_i \left( \; \frac{ y_{i+1} \;-\; y_{i-1} }{ 2 \; \delta x } \; \right) \;+\; Q_i y_i = f_i \\[10pt] \begin{align} y_{i+1} \;-\; 2y_i \;+\; y_{i-1} \;+\; \frac{ \delta x }{2} P_i y_{i+1} \; \;-\; \frac{ \delta x }{2} P_i y_{i-1} \; \;+\; \delta x^2 Q_i y_i \;&=\; \delta x^2 f_i \\[8pt] y_{i+1} \left( 1 \;+\; \frac{ \delta x }{2} P_i \right) \;+\; y_i \left( \;-2 \;+\; \delta x^2 Q_i y_i \right) \;+\; y_{i-1} \left( 1 \;-\; \frac{ \delta x }{2} P_i \right) \;&=\; \delta x^2 f_i \tag{3} \\[15pt] \end{align} $$$$ for \quad i = 1,2,...,n-1 $$Write the coefficients as,
$$ A_i = \left( 1 + \frac{ \delta x }{2} P_i \right) \\[12pt] B_i = \left( -2 + \delta x^2 Q_i \right) \\[8pt] C_i = \left( 1 - \frac{ \delta x }{2} P_i \right) \\[8pt] $$Equation 3 could be represented in matrix form of linear systems
From the given boundary conditions of $$ y_0 \;=\; y\;(x_0) \;=\; y\;(a) \;=\; \alpha \\ y_n \;=\; y\;(x_n) \;=\; y\;(b) \;=\; \beta \\ $$
Let,
$$ x \; = \; \begin{bmatrix} y\;(x_0) \\ y\;(x_1) \\ y\;(x_2) \\ . \\ . \\ . \\ y\;(x_{n-1}) \\ y\;(x_n) \end{bmatrix} \; = \; \begin{bmatrix} y_0 \\ y_1 \\ y_2 \\ . \\ . \\ . \\ y_{n-1} \\ y_n \end{bmatrix} \qquad\qquad b \text{ } = \text{ } \delta x^2 \text{ } \begin{bmatrix} f_0 \\ f_1 \\ f_2 \\ . \\ . \\ . \\ f_{n-1} \\ f_n \end{bmatrix} \; = \; \quad \delta x^2 \begin{bmatrix} \alpha \; \delta x^{-2} \\ f_1 \\ f_2 \\ . \\ . \\ . \\ f_{n-1} \\ \beta \; \delta x^{-2} \end{bmatrix} \\ $$The model expressed as $ A \boldsymbol{x} = \boldsymbol{b} $ for any arbitrary value of $n$ with the boundary conditions in the first and last row,
$$ \begin{align} \\ \begin{bmatrix} 1 & 0 & 0 & 0 & 0 & . & . & 0 & 0 & 0 & 0 \\ C_1 & B_1 & A_1 & . & . & . & . & . & . & . & . \\ 0 & C_1 & B_1 & A_1 & . & . & . & . & . & . & . \\ . & 0 & C_1 & B_1 & A_1 & . & . & . & . & . & . \\ . & . & 0 & 0 & 0 & . & . & . & . & . & . \\ . & . & . & . & . & . & . & . & . & . & . \\ . & . & . & . & . & . & . & C_{n-2} & B_{n-2} & A_{n-2} & 0 & \\ . & . & . & . & . & . & . & 0 & C_{n-1} & B_{n-1} & A_{n-1} \\ . & . & . & . & . & . & . & . & 0 & 0 & 1 \\ \end{bmatrix} \quad \begin{bmatrix} y_0 \\ y_1 \\ y_2 \\ . \\ . \\ . \\ . \\ y_{n-1} \\ y_n \end{bmatrix} \quad &= \quad \delta x^2 \begin{bmatrix} \alpha \; \delta x^{-2} \\ f_1 \\ f_2 \\ . \\ . \\ . \\ . \\ f_{n-1} \\ \beta \; \delta x^{-2} \end{bmatrix} \\ \text{ } \\ A \qquad\qquad\qquad\qquad\qquad \boldsymbol{x} \qquad &= \qquad \boldsymbol{b} \end{align} $$This could be solved as a matrix inversion problem for $\boldsymbol{x}$ to yield,
$$ \boldsymbol{x} \;=\; A^{-1} \boldsymbol{b} $$Given the differential equation,
$$ \frac{d^2y}{dx^2} + 3 \frac{dy}{dx} + 2y = 4x^2 \\[15pt] \text{subject to the boundary conditions,} \quad y\;(1) = 1 \quad \text{and} \quad y\;(2) = 6 \\ $$Comparing to the model
$$ \frac{d^2y}{dx^2} + P\;(x) \frac{dy}{dx} + Q\,(x)\;y = f\,(x) \\[15pt] \text{subject to the boundary conditions,} \quad y\;(a) = \alpha \quad \text{and} \quad y\;(b) = \beta \\ $$Obtain values as,
$$ \begin{align} P(x) &= 3 \\ Q(x) &= 2 \\ f(x) &= 4x^2 \\ \end{align} $$$$ a = 1 \quad\Rightarrow\quad y_0 = \alpha = 1 \\ b = 2 \quad\Rightarrow\quad y_n = \beta = 6 \\ $$Write the coefficients as,
$$ A_i \;=\; 1 + \frac{ \delta x }{2} P_i \quad =\; 1 + \frac{ \delta x }{2} \;(3) \\[12pt] B_i \;=\; -2 + \delta x^2 Q_i \quad =\; -2 + \delta x^2 \; (2) \\[8pt] C_i \;=\; 1 - \frac{ \delta x }{2} P_i \quad =\; 1 - \frac{ \delta x }{2} \; (3) \\[8pt] $$The matrix equation becomes,
$$ \boldsymbol{x} \; = \; \begin{bmatrix} y\;(x_0) \\ y\;(x_1) \\ y\;(x_2) \\ . \\ . \\ . \\ y\;(x_{n-1}) \\ y\;(x_n) \end{bmatrix} \; = \; \begin{bmatrix} y_0 \\ y_1 \\ y_2 \\ . \\ . \\ . \\ y_{n-1} \\ y_n \end{bmatrix} $$Write functions to enable re-use and speed up iterations for multiple values of $n$.
def Aterm( partition_size: float) -> float:
return 1 + (partition_size/2) * 3
def Bterm( partition_size: float) -> float:
return -2 + (partition_size**2) * 2
def Cterm( partition_size: float) -> float:
return 1 - (partition_size/2) * 3
def fx( x_val: float) -> float:
return 4 * (x_val ** 2)
def solve_diff( n: int, lower: int, upper: int, y0: float, y_inf: float,
Afunc, Bfunc, Cfunc, funcf, debug: bool=False) -> np.array:
dx = (upper - lower) / (n+1)
arr_x = np.linspace( lower, upper, n+1)
arr_f = np.hstack( (y0, (dx**2) * np.vectorize(funcf)(arr_x[1:n]), y_inf))
arr_k = (dx**2) * np.vectorize(funcf)(arr_x[1:n])
arr_A = np.zeros( (n+1, n+1))
# assign boundary conditions
arr_A[0,0] = 1
arr_A[n,n] = 1
# assign A, B, C terms
arr_terms = np.array([Cfunc(dx), Bfunc(dx), Afunc(dx)])
# np.array([Afunc(dx), Bfunc(dx), Cfunc(dx)])
for k in range(1,n):
arr_A[k, k-1:k+2] = arr_terms
if debug:
print(f"{dx=}\t{dx**2=}")
print(f"{arr_x=}\n")
print(f"{arr_f=}\n")
print(f"{arr_k=}\n")
print(f"{arr_A.shape=}\n")
print(f"\narr_A=\n{arr_A}\n")
arr_out = np.hstack(( arr_x.reshape(-1,1),
np.linalg.solve( arr_A, arr_f).reshape(-1,1) ))
return arr_out
def to_df( arr: np.array) -> pd.DataFrame:
return pd.DataFrame( arr, columns=["x", "y"])
def style_table(df: pd.DataFrame, columns: int=3, spacing: int=20):
space = "\xa0" * spacing
k, m = divmod(df.shape[0], columns)
output = ""
for j in range(columns):
df_table = df[j*k+min(j, m):(j+1)*k+min(j+1, m)]
df_styler = df_table.style.set_table_attributes("style='display:inline'")
output += df_styler._repr_html_() + space
return output
df10 = to_df(solve_diff( 10, 1, 2, 1., 6., Aterm, Bterm, Cterm, fx))
display_html(style_table( df10, 1) , raw=True)
| x | y | |
|---|---|---|
| 0 | 1.000000 | 1.000000 |
| 1 | 1.100000 | 2.299256 |
| 2 | 1.200000 | 3.288446 |
| 3 | 1.300000 | 4.034290 |
| 4 | 1.400000 | 4.591614 |
| 5 | 1.500000 | 5.005412 |
| 6 | 1.600000 | 5.312546 |
| 7 | 1.700000 | 5.543168 |
| 8 | 1.800000 | 5.721885 |
| 9 | 1.900000 | 5.868738 |
| 10 | 2.000000 | 6.000000 |
df50 = to_df( solve_diff( 50, 1, 2, 1, 6, Aterm, Bterm, Cterm, fx))
display_html(style_table( df50) , raw=True)
| x | y | |
|---|---|---|
| 0 | 1.000000 | 1.000000 |
| 1 | 1.020000 | 1.306593 |
| 2 | 1.040000 | 1.596244 |
| 3 | 1.060000 | 1.869767 |
| 4 | 1.080000 | 2.127942 |
| 5 | 1.100000 | 2.371518 |
| 6 | 1.120000 | 2.601211 |
| 7 | 1.140000 | 2.817710 |
| 8 | 1.160000 | 3.021674 |
| 9 | 1.180000 | 3.213736 |
| 10 | 1.200000 | 3.394503 |
| 11 | 1.220000 | 3.564556 |
| 12 | 1.240000 | 3.724452 |
| 13 | 1.260000 | 3.874727 |
| 14 | 1.280000 | 4.015892 |
| 15 | 1.300000 | 4.148439 |
| 16 | 1.320000 | 4.272837 |
| x | y | |
|---|---|---|
| 17 | 1.340000 | 4.389538 |
| 18 | 1.360000 | 4.498974 |
| 19 | 1.380000 | 4.601560 |
| 20 | 1.400000 | 4.697691 |
| 21 | 1.420000 | 4.787748 |
| 22 | 1.440000 | 4.872095 |
| 23 | 1.460000 | 4.951081 |
| 24 | 1.480000 | 5.025039 |
| 25 | 1.500000 | 5.094290 |
| 26 | 1.520000 | 5.159140 |
| 27 | 1.540000 | 5.219882 |
| 28 | 1.560000 | 5.276797 |
| 29 | 1.580000 | 5.330154 |
| 30 | 1.600000 | 5.380209 |
| 31 | 1.620000 | 5.427211 |
| 32 | 1.640000 | 5.471393 |
| 33 | 1.660000 | 5.512981 |
| x | y | |
|---|---|---|
| 34 | 1.680000 | 5.552192 |
| 35 | 1.700000 | 5.589231 |
| 36 | 1.720000 | 5.624296 |
| 37 | 1.740000 | 5.657576 |
| 38 | 1.760000 | 5.689252 |
| 39 | 1.780000 | 5.719495 |
| 40 | 1.800000 | 5.748471 |
| 41 | 1.820000 | 5.776338 |
| 42 | 1.840000 | 5.803246 |
| 43 | 1.860000 | 5.829340 |
| 44 | 1.880000 | 5.854756 |
| 45 | 1.900000 | 5.879627 |
| 46 | 1.920000 | 5.904078 |
| 47 | 1.940000 | 5.928229 |
| 48 | 1.960000 | 5.952195 |
| 49 | 1.980000 | 5.976083 |
| 50 | 2.000000 | 6.000000 |
For $n=100$,
df100 = to_df( solve_diff( 100, 1, 2, 1, 6, Aterm, Bterm, Cterm, fx))
display_html( style_table( df100, 4, 5), raw=True)
| x | y | |
|---|---|---|
| 0 | 1.000000 | 1.000000 |
| 1 | 1.010000 | 1.156879 |
| 2 | 1.020000 | 1.309337 |
| 3 | 1.030000 | 1.457482 |
| 4 | 1.040000 | 1.601420 |
| 5 | 1.050000 | 1.741253 |
| 6 | 1.060000 | 1.877083 |
| 7 | 1.070000 | 2.009009 |
| 8 | 1.080000 | 2.137128 |
| 9 | 1.090000 | 2.261535 |
| 10 | 1.100000 | 2.382323 |
| 11 | 1.110000 | 2.499583 |
| 12 | 1.120000 | 2.613404 |
| 13 | 1.130000 | 2.723873 |
| 14 | 1.140000 | 2.831077 |
| 15 | 1.150000 | 2.935098 |
| 16 | 1.160000 | 3.036018 |
| 17 | 1.170000 | 3.133918 |
| 18 | 1.180000 | 3.228876 |
| 19 | 1.190000 | 3.320969 |
| 20 | 1.200000 | 3.410273 |
| 21 | 1.210000 | 3.496860 |
| 22 | 1.220000 | 3.580803 |
| 23 | 1.230000 | 3.662172 |
| 24 | 1.240000 | 3.741037 |
| 25 | 1.250000 | 3.817465 |
| x | y | |
|---|---|---|
| 26 | 1.260000 | 3.891523 |
| 27 | 1.270000 | 3.963274 |
| 28 | 1.280000 | 4.032783 |
| 29 | 1.290000 | 4.100111 |
| 30 | 1.300000 | 4.165320 |
| 31 | 1.310000 | 4.228469 |
| 32 | 1.320000 | 4.289615 |
| 33 | 1.330000 | 4.348816 |
| 34 | 1.340000 | 4.406128 |
| 35 | 1.350000 | 4.461605 |
| 36 | 1.360000 | 4.515301 |
| 37 | 1.370000 | 4.567267 |
| 38 | 1.380000 | 4.617555 |
| 39 | 1.390000 | 4.666215 |
| 40 | 1.400000 | 4.713296 |
| 41 | 1.410000 | 4.758846 |
| 42 | 1.420000 | 4.802911 |
| 43 | 1.430000 | 4.845538 |
| 44 | 1.440000 | 4.886771 |
| 45 | 1.450000 | 4.926655 |
| 46 | 1.460000 | 4.965232 |
| 47 | 1.470000 | 5.002544 |
| 48 | 1.480000 | 5.038632 |
| 49 | 1.490000 | 5.073537 |
| 50 | 1.500000 | 5.107299 |
| x | y | |
|---|---|---|
| 51 | 1.510000 | 5.139954 |
| 52 | 1.520000 | 5.171542 |
| 53 | 1.530000 | 5.202099 |
| 54 | 1.540000 | 5.231662 |
| 55 | 1.550000 | 5.260264 |
| 56 | 1.560000 | 5.287942 |
| 57 | 1.570000 | 5.314728 |
| 58 | 1.580000 | 5.340656 |
| 59 | 1.590000 | 5.365757 |
| 60 | 1.600000 | 5.390065 |
| 61 | 1.610000 | 5.413608 |
| 62 | 1.620000 | 5.436418 |
| 63 | 1.630000 | 5.458525 |
| 64 | 1.640000 | 5.479956 |
| 65 | 1.650000 | 5.500741 |
| 66 | 1.660000 | 5.520906 |
| 67 | 1.670000 | 5.540480 |
| 68 | 1.680000 | 5.559488 |
| 69 | 1.690000 | 5.577956 |
| 70 | 1.700000 | 5.595909 |
| 71 | 1.710000 | 5.613372 |
| 72 | 1.720000 | 5.630370 |
| 73 | 1.730000 | 5.646926 |
| 74 | 1.740000 | 5.663062 |
| 75 | 1.750000 | 5.678802 |
| x | y | |
|---|---|---|
| 76 | 1.760000 | 5.694167 |
| 77 | 1.770000 | 5.709180 |
| 78 | 1.780000 | 5.723861 |
| 79 | 1.790000 | 5.738230 |
| 80 | 1.800000 | 5.752308 |
| 81 | 1.810000 | 5.766115 |
| 82 | 1.820000 | 5.779670 |
| 83 | 1.830000 | 5.792991 |
| 84 | 1.840000 | 5.806097 |
| 85 | 1.850000 | 5.819006 |
| 86 | 1.860000 | 5.831735 |
| 87 | 1.870000 | 5.844302 |
| 88 | 1.880000 | 5.856723 |
| 89 | 1.890000 | 5.869015 |
| 90 | 1.900000 | 5.881193 |
| 91 | 1.910000 | 5.893273 |
| 92 | 1.920000 | 5.905271 |
| 93 | 1.930000 | 5.917202 |
| 94 | 1.940000 | 5.929079 |
| 95 | 1.950000 | 5.940917 |
| 96 | 1.960000 | 5.952730 |
| 97 | 1.970000 | 5.964532 |
| 98 | 1.980000 | 5.976336 |
| 99 | 1.990000 | 5.988154 |
| 100 | 2.000000 | 6.000000 |
Obtain the closed form solution of the ODE using SymPy.
y = sp.symbols("y", cls=sp.Function)
x = sp.symbols("x")
eqs = y(x).diff(x).diff(x) + 3*y(x).diff(x) + 2*y(x) - 4*x**2
# boundary condition
bound_cond = {y(1):1, y(2):6}
general_soln = sp.dsolve(eqs, ics=bound_cond)
general_soln
Generate dataframe for plotting.
xval = np.linspace(1,2,600)
y_func = sp.lambdify( x, general_soln.rhs)
df_closedform = pd.DataFrame( np.hstack((xval.reshape(-1,1), y_func(xval).reshape(-1,1) )),
columns=['x','y'])
df_closedform
| x | y | |
|---|---|---|
| 0 | 1.000000 | 1.000000 |
| 1 | 1.001669 | 1.026755 |
| 2 | 1.003339 | 1.053382 |
| 3 | 1.005008 | 1.079881 |
| 4 | 1.006678 | 1.106253 |
| ... | ... | ... |
| 595 | 1.993322 | 5.992169 |
| 596 | 1.994992 | 5.994125 |
| 597 | 1.996661 | 5.996082 |
| 598 | 1.998331 | 5.998041 |
| 599 | 2.000000 | 6.000000 |
600 rows × 2 columns
Generate numerical solution for varying $n$ and plot. (Values of $n$ chosen for visual aesthetics)
ls_df = []
n_vals = [i for i in range(2,10,2)]
for k in n_vals:
ls_df.append( to_df( solve_diff( k, 1, 2, 1, 6, Aterm, Bterm, Cterm, fx)))
fig12 = px.scatter( title="<b>Figure 12 - Varying size of <i>n</i></b>", width=850, height=600
).update_layout(yaxis_range=[0.9,6.1],
xaxis_range=[0.95,2.05])
for i in range(len(ls_df)):
fig12.add_scatter(mode='lines+markers', name=f"n={n_vals[i]:,}",
x = ls_df[i].x,
y = ls_df[i].y,
hovertemplate = 'x: <b>%{x:.6f}</b><br>y: <b>%{y:,.4f}</b>',
line=dict( width=1.5), marker=dict(symbol='triangle-down', size=10)
)
fig12.add_scatter(mode='lines', name=f"Closed Form",
x = df_closedform.x,
y = df_closedform.y,
hovertemplate = 'x: <b>%{x:.6f}</b><br>y: <b>%{y:,.4f}</b>',
line=dict( width=2.5, color='black', dash='dash')
)
# Show spikes
fig12.update_xaxes(showspikes=True, title='x', dtick=0.1)
fig12.update_yaxes(showspikes=True, title='y')
fig12.show()
As $n$ increases, the numerical solution approaches the exact solution. Essentially, $n$ is the number of straight lines to approximate a curve, the larger the value, the finer and better able it is match the exact solution.
Let error be the difference between numerical and general solution.
$$ \varepsilon = \hat{y} - y \\ $$ls_df = []
n1_vals = [2, 5, 10, 25, 50, 100, 200, 500, 1000, 2000, 5000]
for k in n1_vals:
closed_form = y_func( np.linspace(1,2,k+1))
df_num_soln = to_df( solve_diff( k, 1, 2, 1, 6, Aterm, Bterm, Cterm, fx))
df_num_soln['exact'] = closed_form
df_num_soln['error'] = df_num_soln.y - df_num_soln.exact
df_num_soln['error_percentage'] = 100 * (df_num_soln.error / closed_form[1])
ls_df.append( df_num_soln )
ls_df[0].tail()
| x | y | exact | error | error_percentage | |
|---|---|---|---|---|---|
| 0 | 1.0 | 1.00000 | 1.000000 | 0.000000 | 0.000000 |
| 1 | 1.5 | 4.78125 | 5.120806 | -0.339556 | -6.630906 |
| 2 | 2.0 | 6.00000 | 6.000000 | 0.000000 | 0.000000 |
fig13 = px.scatter( title="<b>Figure 13 - Error percentage for varying size of <i>n</i></b>",
width=750, height=800)
for k in range(len(ls_df)):
fig13.add_scatter(mode='lines', name=f"n={n1_vals[k]:,}",
x = ls_df[k].x,
y = ls_df[k].error_percentage,
hovertemplate = 'x: <b>%{x:.6f}</b><br>Error (%): <b>%{y:.5f}</b>',
line=dict( width=1.5))
# Show spikes
fig13.update_xaxes(showspikes=True, title='x', dtick=0.1)
fig13.update_yaxes(showspikes=True, title='Error (%)', tickformat='0.1f')
fig13.show()
Based on the chart above, at around 200 partitions of $x$ the maximum error is less than 1%. And at 2000 partitions, less than 0.1% error. Note that $x=1$ and $x=2$ are the boundary conditions and are fixed.
Negative indicates that the numerical estimate is underestimating the solution.
Solve the following using Monte Carlo integration
$$ \int_{1}^{3} x^2 \,dx \tag{Integral 1} \\[10pt] $$$$ \int_{0}^{\infty} e^{-x^2} \,dx \tag{Integral 2} \\[10pt] $$$$ \frac {1}{\sqrt{2 \pi}} \int_{-\infty}^{\infty} x^4 e^{-x^2/2} \,dx \tag{Integral 3} \\[10pt] $$Monte Carlo method is centered on evaluating definite integrals as expectations.
For example, to estimate $\theta$ given that $$ \theta = \mathbb{E}[f(U)] \quad \text{where} \quad U \sim U(0,1) $$
The expectation can be expressed as an integral of the form $$ \mathbb{E}[f(U)] = \int_{0}^{1} f(u) p(u) \,du $$
where $p(u)$ is the probability density function (PDF) for the uniform distribution $U(,0,1)$ given by
$$ p(u) = \begin{cases} 1 & 0 \le u \leq1 \\ 0 & \mbox{otherwise} \end{cases} $$hence $$ \mathbb{E}[f(U)] = \int_{0}^{1} f(u) \,du $$
There are two fundamental cases depending on the limits of integration.
Case 1: function of $f$ $\text{s.t.} \quad f:[a,b] \rightarrow \mathbb{R} \quad \text{ where } \quad -\infty \lt a \lt b \lt \infty$
For a given integral, $$ I = \int_{a}^{b} f(x)\,dx \quad \\[10pt] $$
transform the domain $[a,b]$ to $[0,1]$ by applying the subsititution,
$$ y = \frac{x-a}{b-a} \quad \Rightarrow \quad x = y \text{ } (b-a) + a $$which gives
$$ dx = dy \text{ } (b-a) $$So now,
$$ \begin{align} I \text{ } &= \text{ } (b-a) \int_{0}^{1} f \left( \text{ } y(b-a) + a \text{ } \right) \,dy \\[8pt] &= \text{ } (b-a) \text{ }\text{ } \mathbb{E} \left[ \text{ } f( \text{ } U (b-a) + a \text{ } ) \text{ } \right] \quad \text{where} \quad U \sim U(0,1) \end{align} $$Case 2: function of $g$ $\text{s.t.} \quad g:[a,\infty) \rightarrow \mathbb{R} \text{ where } \quad -\infty \lt a \lt b \lt \infty$
For a given integral $I$, provided $I \lt \infty$ $$ I = \int_{a}^{b} g(x)\,dx \\[10pt] $$
transform the domain $[a,b]$ to $[0,1]$ by applying the subsititution,
$$ y = \frac{1}{1+x} \quad \Rightarrow \quad x = \frac{1}{y} - 1 $$which gives
$$ dx = \frac{dy} {y^2} $$So now,
$$ \begin{align} I &= \int_{0}^{1} \frac {g \left( \text{ } \frac{1}{y} - 1 \text{ } \right)} {y^2} \,dy \\[10pt] &= \mathbb{E} \left[ \frac{ g \left( -1 + \frac{1}{U} \right) } { U^2} \right] \quad \text{where} \quad U \sim U(0,1) \end{align} $$In all proceeding cases, define error to be the difference between estimate and exact. $$ \varepsilon = \hat{I_1} - I_1 $$
Given the integral,
$$ I_1 = \int_{1}^{3} x^2 \,dx $$the exact solution is,
$$ \begin{align} I_1 &= \left[\frac{x^3}{3}\right]_1^3 \;=\; \left[\frac{3^3}{3}\right] - \left[\frac{1^3}{3}\right] \\[7pt] &= \left[\frac{27}{3}\right] - \left[\frac{1}{3}\right] \\[7pt] &= \frac{26}{3} \end{align} $$Let,
$$
\begin{align}
I_1 &= \int_{1}^{3} f(x) \,dx
\quad \text{where}
\quad f(x) = x^2 \\[7pt]
\end{align}
$$
applying Case 1 to transform from $[1,3]$ to domain $[0,1]$ by the following substitution,
$$ y = \frac{x-1}{3-1} = \frac{x-1}{2} \quad \Rightarrow \quad x = 2y + 1 $$when
$ x = 1$, $ y = 0$
$ x = 3$, $ y = 1$
differentiating $$ dx = 2 \text{ } dy $$
This gives,
$$ \begin{align} I_1 \text{ } &= \text{ } \int_{0}^{1} f( \text{ } 2 \text{ } y + 1 \text{ } ) \, \text{ } 2 dy \;=\; \text{ } 2 \text{ } \int_{0}^{1} f( \text{ } 2 \text{ } y + 1 \text{ } ) \,dy \\[5pt] &= \text{ } 2 \text{ }\text{ } \mathbb{E} \left[ \text{ } f( \text{ } 2U+1 \text{ } ) \text{ } \right] \quad \text{where} \quad U \sim U(0,1) \end{align} $$Next,
exact_I1 = 26 / 3
iters_arr1 = 10**7
def solve_I1(exact: float, iter_arr1: int) -> pd.DataFrame:
np.random.seed(20230414)
arr1 = np.random.uniform(0, 1, iter_arr1)
arr1 = np.vstack((arr1, (arr1*2+1)**2)) # 1 - term
arr1 = np.vstack((arr1, (np.cumsum(arr1[1,:]) / np.arange(1,arr1.shape[1]+1)) )) # 2 - average
# extract values at interval - save RAM and swap
arr1_interval = np.append(1, np.arange(1_000,arr1.shape[1], 1_000))
arr1 = arr1[:,arr1_interval]
arr1 = np.vstack((arr1, arr1[2,:] * 2)) # 3 - average*2
arr1 = np.vstack((arr1, arr1[3,:] - exact)) # 4 - error
arr1 = np.vstack((arr1, (arr1[4,:] / exact) * 100 )) # 5 - error percentage
arr1 = np.transpose(np.vstack((arr1, arr1_interval)))
df = pd.DataFrame( arr1,
columns=['rand','term','average', 'estimate', 'error', 'error_perc','iter'])
return df
df1 = solve_I1(exact_I1, iters_arr1)
# df1.shape
print(f"The closed form solution of Integal 1 has the value {exact_I1:0.6f}\n")
The closed form solution of Integal 1 has the value 8.666667
df1.iloc[:,3:].tail()
| estimate | error | error_perc | iter | |
|---|---|---|---|---|
| 9995 | 8.666773 | 0.000107 | 0.001230 | 9995000.0 |
| 9996 | 8.666760 | 0.000093 | 0.001078 | 9996000.0 |
| 9997 | 8.666764 | 0.000097 | 0.001125 | 9997000.0 |
| 9998 | 8.666788 | 0.000121 | 0.001401 | 9998000.0 |
| 9999 | 8.666793 | 0.000127 | 0.001460 | 9999000.0 |
Table above shows part of the calculated values for Integral 1. Compared to the exact solution, the estimate appears to have converged.
arr1_fig1 = px.scatter( df1, x='iter', y='estimate',
labels={'estimate': 'Estimate', 'iter': 'Iterations'},
title="<b>Figure 14 - Estimate vs Iteration for Integral 1</b><br><sup>(exact solution in red)</sup>",
width=900, height=500,
).update_traces(mode='lines', line = dict(color='blue'),
hovertemplate='Error: %{y}<br>Iteration: %{x}'
).update_layout(yaxis_range=[8.64,8.68])
arr1_fig1.add_hline(y=exact_I1, line_width=1.5, line_dash="solid", line_color="red", opacity=0.7, name="exact")
# Show spikes
arr1_fig1.update_xaxes(showspikes=True)
arr1_fig1.update_yaxes(showspikes=True, tickformat='.2f')
arr1_fig1.show()
Chart above shows the convergence of the Monte Carlo integration for Integral 1. As can be observed, the numerical integration converges relatively fast to some small error after some initial fluctuations. The estimate never seems to truly converge and appears to oscillate around the closed form value. The next chart will show the error percentage to have better sense of performance.
arr1_fig3 = px.scatter( df1, x='iter', y='error_perc',
labels={'error_perc': 'Error (%)', 'iter': 'Iterations'},
title="<b>Figure 15 - Iteration vs Error (%) for Integral 1</b>",
width=900, height=600,
).update_traces(mode='lines', line = dict(color='red'),
hovertemplate='Error (%): %{y:.6%}<br>Iteration: %{x}'
).update_layout(yaxis_range=[-0.10,0.06])
arr1_fig3.add_hline(y=0, line_width=1, line_dash="solid", line_color="black", opacity=0.3)
# Show spikes
arr1_fig3.update_xaxes(showspikes=True)
arr1_fig3.update_yaxes(showspikes=True, tickformat='.2%', dtick = 0.01)
arr1_fig3.update_layout(showlegend=True)
arr1_fig3.show()
The plot above shows the error as percentage against iterations for the Monte Carlo integration. After approximately 7 million iterations, the error is within 1% range. Depending on use case, this may or may not be sufficient accuracy.
Given, $$ \begin{align} I_2 &= \int_{0}^{\infty} e^{-x^2} \,dx \tag{Integral 2} \\[10pt] &= \frac{1}{2} \text{ } \underbrace{ \int_{-\infty}^{\infty} e^{-x^2} \,dx }_{\text{Gaussian integral}} \end{align}\\ $$ since it is known that the Gaussian interval is symmetric about $x=0$.
Notice the Gaussian integral which the answer is known
$$ {\int_{-\infty}^{\infty} e^{-x^2} \,dx} = \sqrt{\pi} $$Therefore $$ I_2 = \frac{\sqrt{\pi}} {2} $$
Let,
$$ I_2 = \int_{0}^{\infty} g(x) \,dx \quad \text{where} \quad g(x) = e^{-x^2} \tag{Integral 2} \\[10pt] $$applying Case 2 to transform from $[0,\infty]$ to domain $[0,1]$ by the following substitution,
$$ y = \frac{1}{1 + x} \quad \Rightarrow \quad x = \frac{1}{y} - 1 $$when
$ x = 0$, $ y = 1$
$ x = \infty$, $ y = 0$
and differentiating gives
$$ dx = - \text{ } \frac{dy} {y^2} $$This gives,
$$ \begin{align} I_1 \text{ } &= - \int_{1}^{0} g( \frac{1}{y} - 1 ) \text{ } \, \frac{ dy }{y^2} \\[7pt] &= \int_{0}^{1} \frac{ g( \frac{1}{y} - 1 )}{y^2} \text{ } \, dy \\[7pt] &= \text{ } \mathbb{E} \left[ \frac{ g( \frac{1}{U}-1) }{U^2} \right] \quad \text{where} \quad U \sim U(0,1) \end{align} $$Next,
exact_i2 = np.sqrt(np.pi) / 2
iter_arr2 = 10**7
def solve_I2( exact2: float, iters: int) -> pd.DataFrame:
np.random.seed(20230414)
arr2 = np.random.uniform(0, 1, iters)
arr2 = np.vstack( (arr2, arr2**2)) # 1 - denominator
arr2 = np.vstack( (arr2, np.exp( -((1/arr2[0,:]) - 1)**2) )) # 2 - numerator
arr2 = np.vstack( (arr2, arr2[2,:] / arr2[1,:])) # 3 - numerator / denominator
arr2 = np.vstack( (arr2, np.cumsum(arr2[3,:]) / np.arange(1,arr2.shape[1]+1) )) # 4- estimate (averaged)
# extract values at interval
arr2_interval = np.append(1, np.arange(2_000,arr2.shape[1], 2_000))
arr2 = arr2[:,arr2_interval]
arr2 = np.vstack( (arr2, arr2[4,:] - exact2)) # 5 - error
arr2 = np.vstack( (arr2, (arr2[5,:] / exact2) * 100)) # 6 - error percentage
arr2_plot = np.transpose( np.vstack( (arr2, arr2_interval )))
df = pd.DataFrame( arr2_plot,
columns=['randuniform', 'denom', 'numer', 'term', 'estimate','error',
'error_percentage','iter'])
return df
df2 = solve_I2(exact_i2, iter_arr2)
print(f"The closed form solution of Integal 2 has the value {exact_i2:0.6f}\n")
The closed form solution of Integal 2 has the value 0.886227
df2.iloc[:,4:].tail()
| estimate | error | error_percentage | iter | |
|---|---|---|---|---|
| 4995 | 0.886315 | 0.000088 | 0.009919 | 9990000.0 |
| 4996 | 0.886318 | 0.000091 | 0.010221 | 9992000.0 |
| 4997 | 0.886318 | 0.000091 | 0.010234 | 9994000.0 |
| 4998 | 0.886317 | 0.000090 | 0.010157 | 9996000.0 |
| 4999 | 0.886318 | 0.000091 | 0.010297 | 9998000.0 |
Table above shows part of the calculated values for Integral 2. Compared to the exact value, the estimate appears to have converged.
arr2_fig1 = px.scatter( df2, x='iter', y='estimate',
labels={'estimate': 'Estimate', 'iter': 'Iterations'},
title="<b>Figure 16 - Estimate vs Iteration for Integral 2</b><br><sup>(exact solution in red)</sup>",
width=900, height=500,
).update_traces(mode='lines', hovertemplate='Estimate: %{y}<br>Iteration: %{x}'
).update_layout(yaxis_range=[0.882,0.888])
arr2_fig1.add_hline(y=exact_i2, line_width=1, line_dash="solid", line_color="red", opacity=0.9)
# Show spikes
arr2_fig1.update_xaxes(showspikes=True)
arr2_fig1.update_yaxes(showspikes=True, tickformat=".3f")
arr2_fig1.show()
Figure above show the estimate performance for increasing iterations. The estimate appears to converge beyond 6 million iterations with some minor oscillations around the exact solution.
arr2_fig3 = px.scatter( df2, x='iter', y='error_percentage',
labels={'error_percentage': 'Error (%)', 'iter': 'Iterations'},
title="<b>Figure 17 - Error (%) vs Iteration for Integral 2</b>",
width=900, height=500,
).update_traces(mode='lines', line = dict(color='red'),
hovertemplate='Error (%): %{y}<br>Iteration: %{x}'
).update_layout(yaxis_range=[-0.1,0.1])
arr2_fig3.add_hline(y=0, line_width=1, line_dash="solid", line_color="black", opacity=0.3)
arr2_fig3.add_hline(y=0.015, line_width=1.5, line_dash="dash", line_color="green", opacity=1)
arr2_fig3.add_hline(y=-0.015, line_width=1.5, line_dash="dash", line_color="green", opacity=1)
# Show spikes
arr2_fig3.update_xaxes(showspikes=True)
arr2_fig3.update_yaxes(showspikes=True, tickformat='.2%', dtick=0.02)
arr2_fig3.show()
The plot above shows the error between the estimate and the exact solution. It appears to be fluctuating up until 5 million iterations before starting to converge and even after that, it does not seem to stay within $\pm$1.5% (green dashed line) and oscillates. Note that for the Gaussian definite integral, a solution is only possible using multivariable calculus, a numerical answer with approximately 1.5% error after 10 million iterations is perhaps passable depending on the situation. Unfortunately, there are no guarantees that more iterations will finally converge. It will be left for further investigation.
Given the integral, $$ \begin{align} I_3 &= \frac {1}{\sqrt{2 \pi}} \int_{-\infty}^{\infty} x^4 e^{-x^2/2} \,dx \\[5pt] \end{align} $$
Apply integration by parts repeatedly onto the integral,
$$ \begin{align} \int_{-\infty}^{\infty} x^4 e^{-x^2/2} \,dx \quad &= \quad \int_{-\infty}^{\infty} x^3 \text{ } \mbox{ } x \text{ } e^{-x^2/2} \,dx \\[10pt] &= \quad \underbrace{ \left[ - e^{-x^2/2} \text{ } x^3 \right]_{-\infty}^{\infty} }_{\lim_{ x\to\infty} \rightarrow 0 } + 3 \int_{-\infty}^{\infty} x^2 e^{-x^2/2} \,dx \\[10pt] &= 3 \text{ } \left( \text{ } \underbrace{ \left[ -e^{-x^2/2} \text{ } x \right]_{-\infty}^{\infty} }_{\lim_{ x\to\infty} \rightarrow 0 } + \underbrace{ \int_{-\infty}^{\infty} e^{-x^2/2} \,dx }_{{\textstyle\unicode{x24B6}}} \text{ } \right) \\[10pt] &= 3 \sqrt{2 \pi} \\[5pt] \end{align} $$For the integral $\textstyle\unicode{x24B6}$,
$$ I = \int_{-\infty}^{\infty} e^{-x^2/2} \,dx = \int_{-\infty}^{\infty} e^{-y^2/2} \,dy \\[15pt] $$Then,
$$ \begin{align} I^2 &= \int_{-\infty}^\infty \int_{-\infty}^\infty e^{ -(x^2+y^2)/2 } \,dx \,dy \\[5pt] &= \int_{0}^{2\pi} \int_{0}^{\infty} e^{-r^2/2} \text{ } r \,dr \,d\theta \\[5pt] &= 1 \text{ } \int_0^{2\pi} \,d\theta \\[5pt] &= 2\pi \\[10pt] \Rightarrow \quad I &= \sqrt{ 2 \pi} \\[5pt] \end{align}\\ $$The exact solution is $$ \begin{align}\\ I_3 &= \frac {1}{\sqrt{2 \pi}} \text{ } \left( \text{ } 3 \sqrt{2 \pi} \text{ } \right) \\[5pt] &= 3 \\[5pt] \end{align}\\ $$
Let,
$$ \begin{align} I_3 &= \frac {1}{\sqrt{2 \pi}} \int_{-\infty}^{\infty} g(x) \,dx \\[10pt] &= \frac {1}{\sqrt{2 \pi}} \text{ } \left[ \text{ } \underbrace{ \int_{-\infty}^{0} g(x) \,dx }_{{\textstyle\unicode{x2460}}}\text{ } + \text{ } \underbrace{ \int_{0}^{\infty} g(x) \,dx }_{{\textstyle\unicode{x2461}}} \text{ } \right] \\[12pt] & \text{where} \quad g(x) = x^4 e^{-x^2/2} \end{align} $$applying Case 2 to transform from $[0,\infty]$ to domain $[0,1]$ by the following substitution,
for integral $ \textstyle\unicode{x2460} $ applying the transformation yields the new domain
$ x = -\infty $, $ y = 0$
$ x = 0 $, $ y = 1$
and for integral $ \textstyle\unicode{x2461} $
$ x = 0$, $ y = 1$
$ x = \infty$, $ y = 0$
and differentiating gives
$$ dx = - \text{ } \frac{dy} {y^2} $$This gives,
$$ \begin{align} I_3 \text{ } &= \frac {1}{\sqrt{2 \pi}} \text{ } \left( \text{ } \int_{0}^{1} g( \frac{1}{y} - 1 ) \text{ } \, \frac{dy}{y^2} \quad - \quad \int_{1}^{0} g( \frac{1}{y} - 1 ) \text{ } \, \frac{dy}{y^2} \text{ } \right) \\[10pt] &= \frac {1}{\sqrt{2 \pi}} \text{ } \left( \text{ } \int_{0}^{1} \frac{ g( \frac{1}{y} - 1) }{y^2} \text{ } \, dy \quad + \quad \int_{0}^{1} \frac{ g( \frac{1}{y} - 1) }{y^2} \text{ } \, dy \text{ } \right) \\[10pt] &= \frac {2}{\sqrt{2 \pi}} \text{ } \left( \text{ } \int_{0}^{1} \frac{ g( \frac{1}{y} - 1) }{y^2} \text{ } \, dy \text{ } \right) \\[10pt] &= \frac {2}{\sqrt{2 \pi}} \text{ } \mathbb{E} \left[ \frac{ g( \frac{1}{U}-1) }{U^2} \right] \quad \text{where} \quad U \sim U(0,1) \end{align} $$Next,
exact_i3 = 3
iter_arr3 = 10**8
def solve_I3( exact: float, iter: int) -> pd.DataFrame:
np.random.seed(20230414)
arr3 = np.random.uniform(0, 1, iter)
arr3 = np.vstack(( arr3, ((1/arr3) - 1) )) # 1 - input to function
arr3 = np.vstack(( arr3, np.square(arr3[0,:]) )) # 2 - U^2
arr3 = np.vstack(( arr3, arr3[1,:]**4 )) # 3 - x^4 (first term)
arr3 = np.vstack(( arr3, np.exp(-np.square(arr3[1,:]) / 2) )) # 4 - exp (second term)
arr3 = np.vstack(( arr3, arr3[3,:] * arr3[4,:] * (1/ arr3[2,:]) )) # 5 - all
# arr3 = np.vstack(( arr3, np.cumsum(arr3[5,:]) / np.arange(1,arr3.shape[1]+1) )) # 6 - averaged
# arr3 = np.vstack(( arr3, (2/np.sqrt(np.pi*2)) * arr3[6,:] )) # 7 - estimate
arr3 = np.vstack(( arr3, (2/np.sqrt(np.pi*2)) * arr3[5,:] )) # 7 - estimate
arr3 = np.vstack(( arr3, np.cumsum(arr3[6,:]) / np.arange(1,arr3.shape[1]+1) )) # 6 - averaged
# extract values at interval
arr3_interval = np.append(1, np.arange(10_000, arr3.shape[1], 10_000))
arr3 = arr3[:, arr3_interval]
arr3 = np.vstack(( arr3, arr3[7,:] - exact )) # 8 - error
arr3 = np.vstack(( arr3, (arr3[8,:] / exact) * 100 )) # 9 - error percentage
arr3_plot = np.transpose( np.vstack( (arr3, arr3_interval )))
df = pd.DataFrame( arr3_plot,
columns=['randuniform', 'input_arg', 'u_square', 'first_term', 'second_term',
'all_terms','averaged', 'estimate', 'error','error_percentage','iter'])
return df
df3 = solve_I3( exact_i3, iter_arr3)
print(f"The closed form solution of Integal 2 has the value {exact_i3:0.6f}")
The closed form solution of Integal 2 has the value 3.000000
df3.iloc[:,7:].tail()
| estimate | error | error_percentage | iter | |
|---|---|---|---|---|
| 9995 | 2.999926 | -0.000074 | -0.002465 | 99950000.0 |
| 9996 | 2.999933 | -0.000067 | -0.002225 | 99960000.0 |
| 9997 | 2.999932 | -0.000068 | -0.002278 | 99970000.0 |
| 9998 | 2.999933 | -0.000067 | -0.002231 | 99980000.0 |
| 9999 | 2.999930 | -0.000070 | -0.002323 | 99990000.0 |
Table above shows part of the calculated values for Integral 3. Compared to the exact value, the estimate appears to have converged.
arr3_fig1 = px.scatter( df3, x='iter', y='estimate',
labels={'estimate': 'Estimate', 'iter': 'Iterations'},
title="<b>Figure 18 - Estimate vs Iteration for Integral 3</b><br><sup>(exact solution in red)</sup>",
width=900, height=500,
).update_traces(mode='lines', hovertemplate='Estimate: %{y}<br>Iteration: %{x}'
).update_layout(yaxis_range=[2.990,3.003])
arr3_fig1.add_hline(y=exact_i3, line_width=1, line_dash="solid", line_color="red", opacity=0.5)
# Show spikes
arr3_fig1.update_xaxes(showspikes=True)
arr3_fig1.update_yaxes(showspikes=True, dtick=0.002, tickformat=".3f")
arr3_fig1.show()
Figure above show the estimate performance for increasing iterations. The estimate appears to converge beyond 6 million iterations with some minor oscillations around the exact solution.
arr3_fig3 = px.scatter( df3, x='iter', y='error_percentage',
labels={'error_percentage': 'Error (%)', 'iter': 'Iterations'},
title="<b>Figure 19 - Error (%) vs Iteration for Integral 3</b>",
width=900, height=500,
).update_traces(mode='lines', line = dict(color='red'),
hovertemplate='Error (%): %{y}<br>Iteration: %{x}'
).update_layout(yaxis_range=[-0.10,0.05])
arr3_fig3.add_hline(y=0, line_width=1, line_dash="solid", line_color="black", opacity=0.3)
arr3_fig3.add_hline(y=0.015, line_width=1, line_dash="dash", line_color="green", opacity=0.9)
arr3_fig3.add_hline(y=-0.015, line_width=1, line_dash="dash", line_color="green", opacity=0.9)
# Show spikes
arr3_fig3.update_xaxes(showspikes=True)
arr3_fig3.update_yaxes(showspikes=True, tickformat='.1%', dtick=0.02)
arr3_fig3.show()
The plot above shows the error between the estimate and the exact solution. Integral 3 required significantly higher iterations before the error was even within 1.5%. It appears to be fluctuating up until 50 million iterations before starting to converge to within $\pm$1.5% (green dashed line).
Monte Carlo integration is a very useful tool to approximate integrals and with the use of the numpy library, only takes 1-2 seconds to calculate depending on the iterations required and the accuracy desired. It is important to note that all 3 integrals seem to converge to within $\pm$1-2% of the exact solution. However, if this was used as part of a larger calcultions, the error accumulation has to be taken into account. Obviously, there is a speed and accuracy trade off. In the numerical estimations, no analysis has been done on stability of the solution. This should be part of the requirement when developing programs to frequently run as some integrals may never converge or contain asymptotic behavior over a certain range.
Further optimization can definitely be done to speed up the code including parallelization and reducing variable assignments.